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TOPICS  IN  MULTIVARIATE  APPROXIMATION  THEORY 


C.  de  Boor 


1.  Approximation  Theory 

In  this  first  lecture,  I  intend  to  give  an  overview  of  what  is 
understood  by  the  term  "Approximation  Theory”.  This  is  both  a  bow 
toward  the  title  of  these  lectures  and  a  survey  of  the  Icinds  of  things 
you  might  reasonably  expect  to  see  covered  in  these  lectures,  albeit 
with  the  special  accent  of  "multivariatr " ,  but  which  I  will  for  the 
most  part  not  cover  at  all.  In  effect,  this  allows  you  to  locate  within 
the  large  scheme  of  things  the  few  specific  items  I  do  cover. 

Approximation  Theory  is  usually  understood  to  deal  first  and  fore¬ 
most  with  best  approxiaation,  or  b.a.  for  short.  This  is  the  task  of 
finding,  given  an  element  x  of  some  metric  space  X  ,  an  element 
m*  from  some  given  subset  M  of  X  for  which 

dist(x,m*)  »  infmgM  dist(x,m)  *:  dist(x,M)  . 

Such  an  m*  is  called  a  b.a.  to  x  from  M  .  In  symbols: 

m*  gfM(x)  . 

Basic  questions  asked  concern: 

Existence:  |$M(x)|  >  0  ? 

Uniqueness:  |1PM(x)|  <  2  ?  More  generally,  |£M(x)  |  *  ? 

Characterization:  How  would  one  recognize  a  b.a.  (other  than  by 
the  brute  force  approach  of  comparing  it  with  all  candidates)?  This  is 
particularly  important  for  the  next  question. 

Construction. 

A  priori  bounds:  What  can  be  said  about  dist(x,M)  based  on  the 
information  that  x  lies  in  some  set  k  ? 

Details  of  the  answers  depend  strongly  on  the  specifics  of  X  , 
dist  ,  and  M  .  Most  commonly,  X  is  a  norrned  linear  space,  such  as 

C(T)  :*  continuous  functions  on  some  locally  compact  metric 
space  T 

(e.g.,  T  ■  (a,b)  or  T  -  R0  or  whatever)  and  the  metric  is  provided 
by  the  norm  on  X  . 
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Existence  requires  that  M  be  closed.  Beyond  that,  it  is  usually 
a  matter  of  local  compactness)  A  minimising  sequence  (mn)  in  N  is 
picked;  this  means  that  m„  6  M  and 

AIq  lx  -  mnl  •  dist(x,N)  . 

Then  (mn  :  n-1,2,...  >  is  bounded,  hence,  by  local  compactness,  has  a 
limit  point  m  in  M  .  For  this  (or  any  other)  limit  point, 

dist(x,M)  <  ix-ml  <  lim  sup  tx-m  I  -  lim  ix-m  I  »  dist(x,M)  , 

n  n 

therefore  m  €  rM(x)  . 

The  standard  example  for  M  is  a  finite  dimensional  linear 
subspace  of  X  ,  e.g., 

*n  i*  polynomials  of  degree  <  n 

as  a  subspace  of  X  ■  C(a,b]  or  L2(a,b]  .  The  desired  local  compact¬ 
ness  is  obvious  for  such  an  M  .  If  M  is  a  nonlinear  subset,  e.g.,  N 
■  *n/*a  i*  rational  functions  of  degree  n  over  m  ,  the  argument  becomes 
more  sophisticated:  The  convergence  notion  used  is  weakened  sufficient¬ 
ly  to  gain  local  compactness  while  not  losing  the  semicontinuity  of  the 
norm  with  respect  to  this  notion  of  convergence. 

A  real  difficulty  in  multivariate  approximation  is  the  fact  that 
it  becomes  reasonable  to  consider  infinite  dimensional  M  .  E.g., 

M  -  C(a,b]  +  C[c,d]  c.  C( [a,b]x(c,d]) 
provides  a  simple  example  of  the  reasonable  attempt  to  approximate  a 
function  of  many  variables  (in  this  case,  two)  by  composition  of  funct¬ 
ions  with  fewer  variables  (in  this  case,  the  sum  of  two  functions  of 
one  variable),  now  even  existence  is  a  nontrivial  matter. 

Uniqueness  and  characterization  involve  a  ball  game,  of  sorts. 
Imagine  the  closed  ball  Br(x)  of  radius  r  around  x  .  Starting 
with  r  ■  0  ,  let  r  grow  until  r  •  r*  :»  dist(x.M)  .  Then 

£m(x)  -MB  .(x)  . 

For  general  M  and  some  x  ,  this  first  touch  may  well  happen  at  two 
or  more  places.  In  such  a  circumstance,  local  uniqueness  and  character¬ 
ization  of  a  local (ly)  b.a.  become  interesting  questions.  The 
interesting  exception  to  this  general  statement  is  provided  by  a 
convex  M  ,  in  which  case  we  have  the  following  picture: 


Geometric  £act.  J^f  M  is  convex  and  m  G  M  and  r  : “ 
then  m  €  FM(x)  iff  it  is  possible  to  separate  M  and  Br(x)  by  a 
closed  hyperplane  through  m  ,  i.e.,  iff  there  exists  X  6  X  for 
which 

*(Bt(x)]  <  Xm  «  X [M]  (1.1) 

The  first  inequality  in  (1.1)  is  equivalent  to 

MBr<0)J  <  X(m-x)  , 

i.e.,  to 

X  *  0  and  IXIr  <  X(m-x)  (  <  IXllm-xl  «  IXlr  )  , 

hence  to 

X  *  0  and  X(m-x)  «  IXllm-xl  . 

This  last  condition  is  called  variously  (  X  *  0  and)  "  m-x  is  an 
extremal  for  X  ■  or  "X  takes  on  its  norm  on  m-x  "  or  "  X  is 
parallel  to  m-x  I  like  this  last  phrase  best  and  therefore  write 
this  condition 

X  | |  m-x 

to  remind  you  of  the  familiar  picture:  In  a  finite  dimensional  setting, 
X  would  simply  be  the  vector  normal  to  the  (separating)  hyperplane 
and  would  point  in  the  same  direction  as  m-x  does;  in  short,  X 
would  be  parallel  to  m-x  . 

As  a  matter  of  convenience,  one  talks  instead  about  the  error 
x-m  .  This  requires  switching  the  sign  of  X  and  so  gives  the 
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Characterization  Theorem.  Let  M  be  a  convex  subset  of  the  normed 
linear  space  X  ,  let  x  €  XXM  ,  and  m  €  M  .  Then 

(i)  m  S?M(x)  <»-■>  3*11  x-n  s.  t.  X [M]  <  Xm  . 

(ii)  _£f  M  is  a  linear  subspace,  then 

m  €  ^M(x)  <— >  3  *  II  x-m  s.t.  X  XM  (i.e. ,  X[M]  -  {0})  . 

Here  is  a  good  exercise  which  can  be  handled  by  entirely  element¬ 
ary  means:  Let  x  6  X  ,  a  nls,  let  X  €  X*  \  0  and  m  e  ker  X  . 

Then,  m  6  X(x)  <mmm>  *<x-m)  »  iXllx-ml  .  For  it,  you  might  want 

to  prove  first  that 

x  6  X  ,  X  €  X*  «■•>  | Xx |  *  lXldist(x,ker  X)  (1.2) 

which  contains  all  the  customary  error  estimates  of  elementary  numeric¬ 
al  analysis. 

Since  X  ||  x-m  and  XXM  together  imply  that  X  ||  x-m'  for 
all  m1  €  FM(x)  ,  nonuniqueness  in  case  of  a  linear  M  is  tied  to  the 
possibility  of  such  a  X  being  parallel  to  more  than  one  element  (of 
the  same  size),  i.e.,  for  the  hyperplane 

(y  e  X  :  Xy  •  1} 

to  touch  the  unit  ball  6^(0)  at  more  than  one  point.  Since  both  M 
and  ^dist(x,M) ace  convex»  having  two  distinct  points  m  and 
m'  in  £m(x)  implies  that  the  whole  line  segment  [m,m']  between 
m  and  m'  is  in  f*M(x)  .  This  says  that  the  unit  ball  must  contain 
line  segments  in  its  boundary,  which  is  the  same  as  saying  that  the 
norm  is  not  strictly  convex.  Put  positively,  strict  convexity  of  the 
norm  (such  as  the  Lp-norm  for  1  <  p  <  •  )  implies  uniqueness  of  b.a. 
from  a  convex  subset. 

For  more  specific  choices  of  X  and  M  ,  the  characterization 
theorem  can  be  made  more  explicit,  to  the  point  where  it  can  be  used 
for  the  construction  of  a  b.a.  For  example,  if  X  »  C(T)  and  is  an 
n-dimensional  linear  subspace  of  X  and  m  €  M  ,  then 
m  €  fM(x)  <■■■> 

3  r<n  ,(wi),(ti)  s.t.  X  i-  I*+1wilti)  II  x"m  *nd  xxm. 

Here,  [t]  denotes  the  linear  functional  of  point  evaluation  at  t  , 
i.e.,  (t ] i f  | — >  f(t)  .  Behind  this  specialization  of  the  general  char¬ 
acterization  theorem  is  the  result  useful  for  Numerical  Analysis  that 
any  linear  functional  on  an  (n-tl)-dimensional  subspace  of  C(T)  has  a 
norm  preserving  extension  of  the  form  • 
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You  will  recognize  in  this  characterization  the  familiar  statement 
that  the  error  x-m  in  a  b.a.  must  take  on  its  norm  at  points  t^, 

. ..,  tr+i  with  r  <  n  and  such  that,  for  some  weights  w^  with 
wi(x-m)(ti)  >  0,  all  i  ,  one  has  rj^w^tj  1  M  • 

It  is  not  difficult  to  see  that  nonunigueness  is  connected  with 
having  r  <  n  here.  (Also,  the  norm  is  not  strictly  convex,  so  we 
would  expect  nonuniqueness  for  some  x  and  M  . )  Recall  that  an 
n-dimensional  subspace  M  of  C(T)  is  called  a  Haar  space  if,  for  any 

distinct  points  t,,  . ..,  t_  in  T  ,  ([t.j)!?  is  linearly  independent 

A  1  r+1 

over  M  .  For  a  Haar  space  M  ,  having  0  *  A  :*  "i^i^  1  M  im“ 

plies  that  r  >  n  .  Having  in  addition  that  \  \ |  x-m  implies  that 

1  | |  x-m'  for  any  other  b.a.  m*  ,  therefore  m  and  m'  must  agree 

at  the  points  tj,  ...,  tr+1  (assuming  without  loss  that  w^  /  0  ,  all 

i)  ,  and,  using  once  more  that  M  is  Haar,  this  implies  that  m  ■  m'  . 

Conversely,  one  can  show  that,  if  M  is  not  Haar,  then  there  are 

functions  with  many  b.a.'s  from  M  . 

This  equivalence  between  uniqueness  and  the  Haar  property  has 

unhappy  consequences  for  multivariate  approximation,  because  of  the 

following 

Fact  (Mairhuber).  I_f  T  is  not  essentially  just  an  interval,  and 
if  dim  M  >  1  ,  then  M  is  not  Haar, 

The  proof  consists  of  a  bit  of  railroading:  Let  (f^,  ...,  fn)  be 
a  basis  for  M  .  Then  det(fj(t^))  is  a  continuous  function  of  the 
n  points  t^,  ...»  tn  .  If  now  T  contains  a  "Y"  ,  i.e.,  a  "fork"  or 
"switch",  then  one  can  continuously  deform  (t^,  t2,  t3,  ...,  tn) 
into  (t2»  t^,  tj,  ...,  tn)  while  keeping  the  tj/s  distinct: 


This  means  that  the  determinant  has  changed  sign  along  the  way,  hence 
must  have  vanished  for  some  choice  of  n  distinct  points. 

The  resulting  nonuniqueness  of  b.a.  in  C(T)  for  multidimension¬ 
al  T  has  produced  a  great  industry  in  uniform  approximation  by  funct¬ 
ions  of  several  variables  and  much  fun  can  be  had.  I  shall  resist  the 
temptation  to  enter  into  details  now,  because  I  am  not  convinced  that 
best  approximation  is  all  that  practical. 


The  question  o£  a  priori  bounds  or  degree  of  approximation  is 
concerned  with 

dist(K,M)  :«  sup  dist(x,M) 
x€K 

A  typical  example  would  be  s  X  ■  C[0,1],  M  -  *  ,  and  K 

n 

{  x  €  X  :  lx"la  <  1  }  .  Actually,  it  is  not  easy  or  even  useful  to  be 
precise  without  getting  simply  the  number  dist(K,M)  .  The  question  of 
degree  of  approximation  comes  into  its  own  when  one  has  given  a  scale 
(Mh)  or  (Mn)  of  subsets  with  h  — >  0  or  n  — >  -  and  then  con¬ 
siders 

EK(h)  :»  dist(K,  Mh) 

as  a  function  of  h  .  One  proves  direct  or  Jackson-type  theorems: 
x  6  K  -«->  dist(x,  Mh)  -  hr 

and  tries  to  demonstrate  their  sharpness,  if  possible,  by  proving 

inverse  or  Bernstein-type  theorems: 

dist(x,  Mh)  ~  hr  -*»>  x/lxl  6  K  . 

Related  is  the  question:  Given  that  x  €  K  ,  is  the  scale  (Mh)  a 

good  choice  for  approximating  x  ?  What  is  one  to  judge  by? 

Kolmogorov  [1936]  proposed  some  time  ago  that  dist(x,M)/dim  M  is  a 

good  measure.  He  introduced 

d_(  K)  :»  inf  dist<K,M)  »:  the  n-width  of  K  . 
dimM  <  n 

While  it  is  not  easy  to  find  an  optimal  subspace,  i.e.,  a  subspace  at 
which  the  infimum  of  the  definition  is  taken  on,  one  would  at  least 
like  an  asymptotically  optimal  scale  (Mn)  ,  i.e.,  a  scale  with  dim  M,, 
■  n  for  which  dist(K,Mn)  ■  0(dR(K))  . 

Once  effort  enters  considerations  of  approximation  (here  in  the 
form  of  the  dimension  of  H  ,  i.e.,  the  degrees  of  freedom  used  in  the 
approximation),  one  can,  of  course,  ask  more:  Is  it  really  worthwhile 
to  construct  best  approximations,  particularly  when  a  near-best  approx¬ 
imation  is  cheaply  available?  Here  we  call  the  linear  map 
A:X  — >  M  a  near-best  approximation  scheme  if 

3  const  V  x  6  X  ix-Axl  <  const  dist(x,M)  . 

Here,  const  -  1  would  be  best  possible.  Ax  would  then  be  a  b.a.  for 
every  x  .  This  does  happen  in  inner  product  spaces  but  hardly  anywhere 
else  except  in  very  special  circumstances.  In  any  event,  such  A  is 
necessarily  a  linear  projector  (with  ran  A  *  M  )  since  the  inequality 
implies  that  A|M  »  1  .  These  approximation  schemes  are  the  topic  of 
the  second  lecture. 


There  ere  many  books  on  Approximation  Theory  available.  One  of  the 
most  striking  is  Lorentz  [1966}  .  Akhiezer  [1967]  summarizes  the  clas¬ 
sical  part.  Both  Cheney  [1966]  and  Rivlin  [1969]  provide  a  careful 
modern  introduction  to  the  field  while  the  two  volumes  of  Rice  [1964, 
1969]  bring  quite  a  bit  of  additional  material,  especially  on  approxim¬ 
ation  from  a  nonlinear  M  .  Powell  [1981]  and  Schfinhage  [1971]  each 
give  a  very  interesting  view  of  the  subject. 


2.  Linear  Interpolation 


Linear  projectors  arise  from  interpolation,  as  I  intend  to  make 
clear  in  this  lecture.  Fortunately  for  me.  Ward  Cheney  is  with  us  who 
has  spent  a  good  part  of  his  professional  life  studying  linear  project¬ 
ors.  He  will  no  doubt  be  ready  to  answer  all  questions  left  over  after 
(or  raised  by)  this  lecture.  Look  for  his  publications  (e.g.,  Morris  a 
Cheney  (1974]).  See  also  the  excellent  book  by  Davis  (1963]. 

The  setup  is  quite  simple:  We  have  a  linear  space  X  of  functions 
on  some  domain  T  and,  correspondingly,  the  linear  space  X*  of  lin¬ 
ear  functionals  on  X  .  We  have  given  f^  ,  . ..,  fm  €  X  and 

li  *  • ..,  X  €  X'  and  consider  the 
l  n 

Task:  Given  g  6  X  ,  construct  Pg  : *  £™a(j)f^  such  that  Pg 

interpolates  to  g  at  ,  . ..,  XR  ,  i.e., 

X^Pg  -  X^g  ,  i«l,...,n. 

For  example,  the  specifics:  T  *  (a,b]  ,  X  »  C[a,b]  ,  fj  * 

(jJ-1  and  X ^  =  (t^J  describe  the  task  of  polynomial  interpolation. 
Altering  this  to  X^f  | — >  /  f (x)f ^xjdx  ,  all  i  ,  and  choosing  m  *  n 
describes  least-squares  approximation  by  polynomials. 

Our  first  observation  is  that  this  task  does  not  depend  on  the 
individual  functions  f^  ,  fm  nor  on  the  individual  linear  funct¬ 

ionals  X^  ,  . ..,  Xr  ,  but  only  on  their  spans 

F  span  (fj™  :-  (I™  t(j)  i  t  6  rf"} 

and 

A  :■  span  (x^)”  • 

This  is  obvious  for  the  f^'s  since  the  very  task  is  stated  in  terms  of 
their  span.  As  to  A  ,  observe  that 

-  *Ah  ,  i-l,...,n  <--«>  V8  €  Rn  (l@(i)Xi)g  «  (lB(i)Xi)h  . 

We  use  the  abbreviation 

LIP(F, A) 

for  the  Linear  Interpolation  Problem  given  by  F  and  A  ,  i.e.,  for 
the 

Task:  Given  g  6  X  ,  find  Pg  6  F  s.t.  g  ■  Pg  on  A  . 

Here,  F  and  A  are  understood  to  be  linear  subspaces  (finite  dimens¬ 
ional)  of  X  and  X'  ,  respectively.  We  call  LIP(F,A)  correct  if  the 
task  has  exactly  one  solution  for  every  g  6  X  . 
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Now,  having  just  gotten  rid  of  the  f^'s  an<^  *i's  »  Xs  con“ 

venient  to  reintroduce  them,  in  a  possibly  refined  form:  Let  (f^  be 
a  basis  for  F  and  let  (x.^  be  a  basis  for  A  .  Then 
Ij<»{j)fj  solves  LIP  ( F ,  A )  for  given  g  <===> 

a  €  Rm  solves  I^(xif^)a(j)  *  X^  ,  i=l,...,n. 

We  conclude 

Lemma,  (i)  LIP(F,A)  is  correct  <===>  the  Gramian  G 

«-  (Xifj^i!l,j2l  is  inverUble. 

(ii)  LIP ( F , A )  is  correct  ===>  Pg  =  Ej<Mj)fj  with 
a  -  A"1(Aig)  . 

The  proof  is  linear  algebra:  Since  (f^)  is  linearly  independent, 
uniqueness  is  equivalent  to  having  A  1-1.  Since  (X^)  is  linearly 
independent,  existence  is  equivalent  to  having  A  onto.  Note  that 
correctness  implies  m  =  n  . 

The  map  P  defined  by  such  a  correct  LIP  is  linear  (as  a  compos¬ 
ition  of  the  linear  maps  g  1 — >  (x^g)  | — >  a  =  G  1(xig)  ) — >  Ija(j)fj  ). 
Also,  by  uniqueness,  P|F  *  1  ,  hence  P2  *  P  >  showing  that  P  is  a 
linear  projector.  Its  range  is 

ran  P  -  F  «  (x  6  X  :  Px  =  x) 
while  its  kernel  or  nullspace  is 

ker  P  *  (x  €  X  :  Xx  *  0,  all  X  6  A)  *:  Ax  *  ran(l-P)  . 

The  customary  view  of  a  linear  projector  is  that  it  provides  a  direct 
sum  decomposition: 

X  *  Px  +  (l-P)x  . 


Px  is  the  projection  of  x  onto  ran  p  along  ker  P  .  I  prefer  to 
think  of  P  as  given  by  F  =  ran  P  and  its  interpolation  conditions 
A  ■  (v  6  X'  :  uP  =  u)  =  ran  P*  . 

This  stresses  the  fact  that  Px  is  the  unique  element  in  ran  P  which 
agrees  with  x  on  A  . 

The  construction  of  Pg  involves,  off  hand,  the  solution  of  the 
linear  system  Ga  ■*  (x^g)  .  This  can  be  viewed  as  switching  over  to  the 


new  basis  (A|)  for  A  with 

Such  a  basis  is,  by  its  construction,  dual  to  (f^)  ,  i.e., 

xifj  “  sii  • 

Another  possibility  is  the  Lagrange  approach:  Switch  over  to 

fj  *- 

so  that  now  (A^)  “  1*  Vet  another  possibility  is  the  Newton 
approach:  If  possible,  switch  over  to 

fj  :-  I^li.jlfj  ,  A*  :-  IjL"1U,j)Aj 
with  LR  a  triangular  factorization  for  G  ,  giving  again  (AJf^)  *  1 
.  In  this  last  approach,  we  would  need  invertible  for 

r*l,2,...  .  Equivalently,  we  would  need  the  LIP(Fr,Ar>  with  Fr  : »  span 
(f1,...,fr)  ,  Ar  span  (A1,...,Ar>  to  be  correct,  giving  rise  to  a 
projector  Pr  ,  r*l,2,...,n  .  In  these  terms,  suitable  bases  for  F 
and  A  can  be  constructed  bootstrap  fashion: 


Xi(1-pi-l> 


jJi'W  /xi£i 


which  is,  in  effect.  Gauss  elimination  without  pivoting. 

The  example  inspiring  all  this  terminology  is,  of  course,  polynom¬ 
ial  interpolation,  mentioned  earlier,  in  which  fj  -  ()^_1  and 
Xi  “  •  The  Lagrange  approach  leads  to 

£j  lt*“>  ,i#j«t-tj)/(ti-tj) 
while  the  Newton  approach  leads  to 

fj  :t  >  "i^tt-t.)  ,  AJ  -  [tx, . . .  ,ti]  :« 
the  divided  difference  at  tj,  . ..,  t^  ,  and  thence  to  the  Newton  form 


Pg  **  l t^ , . . . , t ^ ) g  H  {•  —  tj)  . 

Another  well  known  example  is  specified  by  :  X  *  C(T),  F 
*  span  (f^)  ,  of  dimension  n  ,  and 

A^:g  )— >  /Tfig  ,  i>l,...,n  . 

Now  P  is  Least-squares  approximation,  and  the  Newton  approach  is,  in 
this  instance,  called  Graa-Schmidt  orthogonalization. 


We  are  interested  in  linear  projectors  because  they  provide  near¬ 
best  linear  approximation  schemes.  Explicitly,  we  have 

Lebesgue's  inequality:  Ig  -  Pgl  <  11-PI  dist(g,F) 
in  case  X  is  a  normed  linear  space  (■:  nls)  and  P  is  bounded.  In 
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fact,  we  have  a  bit  more: 

|ug  -  uPg|  <  dist(u,  A)  ll-Pl  aistt9,  F)  ,  y  €  X*,  g  6  X  .12.1) 

For  the  proof,  note  that,  for  i  6  A  and  f  6  F  ,  A(l-P)  *  0,  (l-P)f 

-  0  ,  so 

|ug-uPg|  «  | <v-i)(l-P) (g-f) |  <  ly-Ai ll-Pl ig-f i  . 

Now  take  the  infimum  over  A  6  A  and  f  6  F  . 

The  inequality  (2.1)  is  important  for  the  rule  makers  who 

customarily  approximate  ug  by  yPg  ever  since  Newton  proposed  this 
for  ug  ■  / g  and  P  polynomial  interpolation.  The  variational 
approach  to  splines,  particularly  important  for  the  understanding  of 
Duchon's  multivariate  'thin  plate'  splines  (Duchon  [1976],  [1977], 
Meinguet  [1979]),  takes  off  from  this  setup.  See  the  epilogue. 

The  basic  inequality  (2.1)  raises  the  two  questions: 

(i)  Is  P  bounded?  (ii)  How  big  is  IPI  (or,  ll-Pi  )? 

The  following  two  lemmas  give  answers  of  sorts. 

Lemma  2.1.  Let  p  be  given  on  the  nls  x  b£  LIP(F,A).  Then 
P  is  bounded  iff  ACx*  (:*  continuous  linear  functionals  on  x  ). 

Proof.  "■=>■  For  all  A  €  A  ,  A  »  AP  »  (A|F)op  with  A  |  F  cont¬ 
inuous  since  dim  F  <  •  . 

*<■■"  We  can  write  P  -  1. *  hence  ipi  <  E“if!iiA|l  <  -  . 

Ill 

Lemma  2.2.  IPI  *  sup  inf  lAlifi/|Af|  . 

f€F  A6A 

Proof.  For  any  linear  map  A  , 

■Al  a  sup  lAxl/lxl  >  sup  sup  lAxl/lx-yl 
x(*kerA  x£kerA  ySkerA 

-  sup  l Axi/dist (x ,kerA)  -  sup  inf  lAxllAl/|Ax|  . 
xgkecA  x0kerA  AlkerA 

If  now  A  is  the  linear  projector  P  ,  then  x  »  Px  +  (l-P)x  and  Px 
€  F  -  ran  p  while  (l-p)x  €  ker  p  .  Further,  A  1  ker  P  iff 
A  6  A  (■  ran  P' )  ,  using  the  fact  that  ran  P  *  {x  :  x  «  Px)  ,  hence 
A  1  ker  P  — *>  A  l  ran(l-P)  =■>  a  *  up  *«>  A  €  ran  p'  =  A  . 
Therefore , 

IPI  »  sup  inf  ipx!IAI/|Ax|  ■  sup  inf  ixllAl/|Ax|  .  ||| 

xjSkerP  AikerP  x€F  ASA 
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3.  The  tensor  product  construct 


Our  first  foray  into  a  multivariate  setup  is  by  tensor  products. 
This  construct  is  of  limited  use.  Yet  when  it  can  be  employed ,  it  is  so 
efficient  that  it  is  worth  some  effort  to  bring  a  given  approximation 
problem  into  this  form,  if  it  can  be  done  at  all.  Somewhat  surprising 
approximation  theoretic  advantages  of  tensor  products  are  discussed  in 
de  Boor  t  Devore  (1981). 

The  mathematics  is  quite  simple,  yet  papers  still  appear  which 
look  forbidding  and  needlessly  complicated  since  they  do  not  make  full 
use  of  the  fact  that,  when  dealing  with  tensor  products,  everything  is 
essentially  univariate,  even  the  computer  programs. 

Naively,  the  tensor  product  of  two  univariate  linear  approximation 
schemes  or  projectors  p  and  Q  is  formed  as  follows.  For  each 
fixed  y  ,  the  linear  projector  P  is  applied  to  the  y-seetion 

hy  :*  h( • ,y) 

of  the  bivariate  function  h  ,  giving  the  univariate  function 

Phy  «  «( j shy ) f j 

in  which  the  coefficients  o(j)  ■  a(j;hy)  depend  on  hy  ,  hence  on  y. 
Then  Q  is  applied  to  each  of  these  coefficient  functions 

Cj»y  \ — >  «( j;hy)  , 

thus  obtaining  their  univariate  approximations 

QCj  -  £k6(Jcsj)gk, 

with  (gk)  a  basis  for  ran  Q  .  Altogether,  this  gives  the  approximat¬ 
ion 

(Bh)(x,y)  «  E  #(ksj)gk(y)f ■ (x) 

to  h  .  3'*  3 

Several  questions  are  immediate:  What  is  the  nature  of  R  ?  is 
R  a  projector?  What  is  the  corresponding  LIP?  If  we  first  use  Q  on 
each  x-section  of  H  and  then  P  on  the  resulting  coefficient  funct¬ 
ions,  would  the  resulting  approximation  be  again  Rh  ? 

We  consider  the  last  question  first,  since  its  answer  supplies 
also  the  answer  to  the  other  questions.  At  its  root  is  the  question  of 
whether  the  two  operators  P91  and  19Q  defined  by 

(P«l)h  s  (x,y>  >  (Phy)(x)  ,  (l«Q)h  :  <x,y)  |— >  (Qhjty) 

commute.  Since  both  p  and  Q  are  linear  projectors,  this  question  is 
settled  once  we  know  that  X«1  and  19u  commute  for  each  X  €  ran  P' 
and  each  v  €  ran  Q'  .  This  is  dealt  with  next  in  careful  (and  perhaps 
boring)  detail. 
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Suppose  that  xr  is  a  Is  of  functions  on  some  domain  Tr  ,  r-1,2. 
The  tensor  product  X1®X2  of  the  two  linear  spaces  XA  and  X2  i» 
customarily  defined  as  the  dual  of  the  linear  space  of  all  bilinear 
functionals  on  X^Xj  .  Since  I  intend  to  use  this  concept  only  in  the 
context  of  function  spaces.  I  prefer  the  following  definition:  For  xr 


e  x. 


X1®X2 


T1*T2 


— >  R 


(W  > 


X1 ^  fci >x2 1  fc2  > 


defines  a  function  on 


Ti*T2 


called  the  tensor  product  of 


Further. 


x1*x2 


span{  x1«x2 


xr  6 


xr) 


VT2 


It  is  easy  but  essential  to  verify  that  the  map 

XjXXj  — >  XX«X2  : 

is  bilinear.  This  implies  that 

{  *  x— 1  € 


^X1,X2^  I — > 


R 


X1®X2 


X1®X2 


V 


£j  —  ’  "rj 

(i.e.,  the  scalars  in  the  linear  combinations  making  up 

absorbed ) . 

Prime  examples  are  the  spaces 


and 


P_®P. 


n 


polynomials  in  two  variables,  of  degree  < 
and  of  degree  <  n  in  the  second. 


m 


Xj8X2  can  be 


in  the  first 


and 

rf“®Rn  *  h”1*11  *  m-by-n  matrices  . 

Correspondingly,  one  defines  the  tensor  product  A«u  of  A  €  Xj^ 
and  11  €  X£  as  a  linear  functional  on  X1®X2  by  the  rule 
A«U  j  XjSXj  — >  R  s  Ij  Jj  (Xuj  K“vj  )  • 

This  is  so  obviously  a  linear  functional  on  X^Xj  that  it  is  easy  to 
overlook  the  only  nontrivial  (though  certainly  not  very  deep)  point 
about  this  definition,  viz.  whether  it  is  a  definition.  The  problem  is 
that  the  rule  for  the  value  of  A«u  at  w  »  I^u^av^  makes  explicit 
use  of  the  particular  representation  of  w  mentioned.  Elements  of 
Xj«X2  have  many  different  representations.  For  example,  if  u  =  a+b 
and  v  ■  2c  ,  then 

u®v  ■  a®v  +  b®v  «  a®v  +  (2b)«c  »  a 8c  +  aSc  +  b«v  =  ... 

We  should  therefore  have,  correspondingly,  that 
(lu)  (yv)  «  (Aa)(pv)  (Ab)(xv)  = 

■  (Aa)(xv)  +  (2Ab)(wc)  * 

■  (Aa)(wc)  +  (Aa)(vc)  +  (xb)(uv)  ■  ... 

We  must  therefore  verify  that  the  number  Ej  (Au..  )  (uv ^  )  depends  only  on 
A  ,  v  ,  and  w  :=  Jjuj®vj  '  and  not  on  the  particular  representation 
for  w  used ,  i.e.,  not  on  the  individual  Uj's  and  Vj's.  For  this. 
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let 

Then 


v  £j(Xuj)vj  . 

V  €  x2  ,  and 

v(t2)  -  (l<Xuj)vj)(t2)  -  I(Xuj)vj(t2) 


■  xClu^v^Ctj) )  -  Xw(*,t2) 

showing  that  v  is  the  X-section  of  w  , 

V  -  wx  . 

In  particular,  v  depends  only  on  w  and  X  .  On  the  other  hand, 
Ej(Xu;j)(uvj)  -  wwA  . 

Note  that  we  get  by  symmetry  that  also  Ij(Xu^)(gVj)  »  Xw^  , 

showing  the  hoped-for  commutativity. 

We  are  ready  to  define  the  tensor  product  of  two  LIPs  (did  the 

T 

Tiny  Tim  craze  ever  catch  on  here?):  Assume  that  Fr  c.  Xr  C.  R  r  and 
ArCXr*  give  rise  to  the  correct  LIP(Fr,Ar)  with  corresponding 
invertible  Gramian 


G  :• 
r 


(XriCrj)i!j-l 


Pick  some  space  W  of  functions  on  T^Tj  which  contains 
Fj®F2  .  Further,  pick  some  6  W*  for  which 

wij|F  *  (Xli*X2j)|F  '  ail  ' 

and  set  A 


F  :« 


span( ) 


Then:  (i)  (fAj)  :*  (f 


li*f2j) 


is  a  basis  for  F  . 

(ii)  The  LIP( F , A )  is  correct. 

(iii)  The  corresponding  interpolant  Pw  to  w  6  W  can  be  comput- 


ed  as 


Pw  ■  E  r (i,j) 
i.i  w 

with 

rw  ■  Gl’1S#(G2’1>T  Lw(i'3)  *■  vijw  »  i»j  *_ 

For  the  proof,  any  w  €  F  can  be  written  w  -  I^Hi,])^ 

some  suitable  coefficient  matrix  r  .  From  this,  we  compute  that 
LwU.s)  vrsw  - 

E  G, (r, i)  r<i, j)  G-(s,?‘  -  -  -  ~  T 

i»j  1  2 


,j)  ■  gx  r  g2 


(3.1) 

for 


-1  T  -1 

This  shows  that  r  ■  G^  LW(G2  )  ,  i.e.,  r  depends  only  on  w  . 

This  proves  (i).  Further,  it  shows  that  (x^w)  »  L  <■■■> 

T  ■  Gj"XL(G2t)“x  ,  thus  proving  (ii)  and  (iii). 
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A  simple  example  is  provided  by  the  tensor  product  of  cubic 
Hermite  interpolation,  i.e.,  cubic  interpolation  at  the  four  points 
a,  b,  b,  with  osculatory  parabolic  interpolation,  i.e.,  parabolic 
interpolation  at  the  three  points  c,  c,  d  .  The  above  description 


a, 


leaves  considerable  freedom  in  the  choice  of  the 


choice  would  be 


wijw 


i\i  1 

1  W<a,c) 

2  wx(a,c 

3  w(b,c) 

4 


wx(b,c) 


Wy  ( a , c ) 

wXy(a,c) 

Wy(b,c) 

WXy{b,0 


«  3 
w(a,d) 
wx(a,d) 
w(b,d) 
wx(b,d) 


A  natural 


and  this  would  require  the  function  space  W  to  consist  of  functions 
defined  (at  least)  on  the  rectangle  [a,bj*[c,d]  (assuming  that  a  < 
b  and  c  <  d  )  and  to  have  first  derivatives  and  the  first  cross  deri¬ 
vative  (in  a  pointwise  sense).  Further,  the  desired  commutativity  would 


require  that  w 


for  all  w  €  w 


"xy  yx 

The  computational  advantage  of  this  construct  over  other  means  of 
approximation  is  considerable.  It  generally  takes  0(N3/3)  operations  to 
solve  the  linear  system  for  the  coefficients  of  the  solution  of  a  LIP 
using  N  degrees  of  freedom.  This  number  can  be  reduced  considerably 
in  a  univariate  setup  (such  as  in  spline  approximation)  through  the  use 
of  special  bases  which  make  the  resulting  system  banded.  This  is  much 
harder  in  a  multivariate  context.  In  any  event,  if  the  LIP(Fr,Af)  in¬ 


volves 

volves 


degrees  of  freedom,  r-1,2  ,  then  their  tensor  product  in- 


Yet,  using  (3.1),  one  can  compute  the 

-2  +  n23/3) 


0(n13/3  +  ni2nj 


+  n.n 


ln2 


ijn2  degrees  of  freedom, 
requisite  nj^  coefficients  in 
operations,  since  it  only  involves  solving  n2  systems  of  order  n^ 
with  the  same  coefficient  matrix,  and  n^  systems  of  order  n2  with 
the  same  coefficient  matrix.  Further  savings  are  possible  because  this 
reduction  of  the  computations  to  the  univariate  context  provides  the 
opportunity  to  make  use  of  whatever  savings  are  available  in  that  con¬ 
text.  if,  for  example,  the  univariate  schemes  give  rise  to  banded  syst¬ 
ems,  then  their  tensor  product  can  be  constructed  in  0(n^n2)  operat¬ 
ions. 


This  leads  to  a  point  made  in  de  Boor  [1979] ,  that  it  is  possible 
to  form  the  "tensor  product"  of  the  computer  programs  for  the  solution 
of  the  univariate  "factor "-problems.  Typically,  one  has  available  pro¬ 
grams  INTERNS, M, A)  which  take  the  input  vector  B  *  ^ri^l*"  w*th 
M  *  nr  ,  and  produce  from  it  the  desired  coefficient  vector  A  *  a  for 
the  interpolant  Prg  >»  E  a(j)fr^  to  g  .  Two  changes  are  required  to 
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make  such  a  program  amenable  to  tensor  product  computations:  One  ex¬ 
tends  it  to  allow  B  to  be  a  collection  of  N  input  vectors,  i.e.,  to 
allow  B  to  be  an  M-by-N  matrix  having  these  N  data  vectors  as  its 
columns;  this  requires  that  also  N  be  supplied  on  input.  Correspond¬ 
ingly,  the  output  A  becomes  a  collection  of  N  coefficient  vectors, 
i.e.,  a  matrix.  But,  for  reasons  that  will  be  obvious  in  a  moment,  the 
resulting  coefficient  vectors  should  form  the  rows  of  that  matrix, 
i.e.,  A  should  be  an  N-by-M  matrix.  This  gives  the  extension 
INTERj. ' (B,M,N,A)  .  With  it,  we  can  solve  the  tensor  product  of  the 
two  LIP(Fr,Ar>  ,  r-1,2  ,  by  the  two  calls 
CALL  INTERj' (Ly.nj.njjA) 

CALL  1NTER2‘ (A,n2,n1,r) 

The  two  programs  could  even  be  the  same  as,  e.g.,  in  bicubic  spline 
interpolation.  Further,  once  such  extended  programs  have  been  written, 
it  is  possible  to  carry  out  tensor  product  interpolation  using  more 
than  two  factors.  Finally,  this  formation  of  program  "tensor  products" 
is  also  helpful  in  the  evaluation  or  other  manipulation  of  the 
interpolant.  In  any  event,  the  detailed  programming  effort  takes  place 
at  the  univariate  level,  just  as  in  the  mathematical  analysis  of  the 
construct. 

The  error  is  easily  obtained  formally.  Writing  Er  :»  1  -  Pr  for 
the  error  in  the  linear  approximation  scheme  Pr  ,  we  find 

1  -  (P^EjWPj+Ej)  -  piep2  +  P1,E2  *  E1®P2  +  E1®E2  *  (3,2) 

This  shows  the  error  operator  for  the  tensor  product  scheme  Pj®p2  to 
be  a  sum  of  the  univariate  errors.  The  order  of  approximation  is  there¬ 
fore  no  better  than  the  worse  of  the  two  univariate  schemes. 

Now  note  that  PrEr  -  ErPr  *  0  .  This  implies  that  any  sum  of 
terms  from  the  right  hand  side  of  (3.2)  gives  a  linear  projector.  In 
particular,  Gordon  (1969) 1>2  has  proposed  the  use  of  the  socalled 
Boolean  sum 

P1®P2  S“  Pl®1  +  1®P2  "  P1®P2  *  P1®P2  +  P1®E2  +  E1*P2  * 

The  resulting  approximation  scheme  is  called  blending  since  it  uses 

interpolation  conditions  of  the  form  1B1  and  lBu  ,  hence,  in  its 
simplest  form,  matches  information  along  certain  lines  parallel  to  the 
axes  and  so  constructs  a  surface  by  "blending"  together  certain  curves. 
For  blending,  the  error  is  the  product  of  the  univariate  errors.  This 
improvement  over  the  tensor  product  is  bought  at  a  high  price:  An  in¬ 
finite  amount  of  information  about  the  function  to  be  approximated  is 
required.  Gordon  has  dealt  successfully  with  this  problem  by  proposing 
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that  one  first  use  a  relatively  dense  but  finite  amount  of  information 
to  construct  good  approximations  to  the  required  curves  and  then  use 
these  approximations  in  the  final  construct. 


4.  Multivariate  polynomial  interpolation 

I  begin  with  a  review  of  the  standard  notation  concerning  polynom¬ 
ials  in  m  variables.  The  notation  is  designed  to  make  it  all  look 
just  as  in  the  univariate  case.  The  general  polynomial  of  total 
degree  <  k  is  ,  by  definition,  any  linear  combination 

x  >  r|«|<k  A«  x“ 
of  the  monomials  ()a  with  |a|  <  k  .  Here, 

x*  x(l)a<1)*  ...  •  x(m)a(^1, 

and  the  length  |a|  of  the  integer  vector  a  is  defined  by 
|a|  j-  lalx  -  a(l)  ♦  ...  +  a(m) 

if,  as  we  assume,  all  the  components  of  a  are  nonnegative.  For  such 
an  index  vector,  one  sets 

al  :»  a( 1) 1 •  ...  *a(m)l 
and  thereby  recovers  the  binomial  formula 

(x  +  y)a  -  I  (®)  xay0  . 

8<a  B 

The  partial  ordering  used  here  is  componentwise: 

8  «  a  :»  for  all  i,  8(i)  <  a(i)  . 

This  gives  Leibniz'  formula 

D*(fg)  -  Z  (®)  (DSf)(Da~Sg) 

8<a  p 

for  the  derivative  of  a  product.  Here, 

D°  *-  D*(m) 

l  m 

with  the  partial  derivative  with  respect  to  the  ith  argument.  More 

generally,  p (D)  is  the  constant  coefficient  differential  operator 

p(D)  :=  Z  AD* 

«  CL 

in  case  p  is  the  polynomial 

P  *  JaV)a  * 

For  the  special  linear  polynomial 

p:x  | — >  x*y  Zix(i)y(i)  , 

we  write 

Dy  Z  y(i)Di 
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instead  of  D*y  for  the  resulting  (unnormalized)  derivative  in  the 
direction  of  y  . 

In  one  variable,  it  is  convenient  to  talk  about  Pk  ,  the  linear 
space  of  polynomials  of  order  k  ,  i.e.,  of  degree  <  k  ,  since  its  dim¬ 
ension  is  k  and  the  optimal  approximation  order  from  Pk  achievable 
on  an  interval  of  length  h  is  h*  .  In  several  variables,  approxim¬ 
ation  order  continues  to  be  linked  to  (total)  polynomial  order,  but  the 
dimension  and  other  interesting  quantities  are  more  easily  expressed  in 
terms  of  (total)  degree  rather  than  order.  For  this  reason,  I  will 
concentrate  on  the  linear  space 

•k  -  V*”) 

of  all  polynomials  of  (total)  degree  <  k  in  m  variables.  It  is  not 
difficult  to  see  that 

dim  -  (".*)  • 

Indeed,  the  rule 

«(r)  »-  i(r)  -  i(r-l)  -  1  ,  r-l,...,m 
sets  vip  a  1-1  correspondence  between 

{«  e  :  |  a  |  <  k) 

i«d  the  set 

(I  £  (l,...,«+k)  :  |I|  ■  m) 

of  iW-rdinality  ("**)  if  we  let  i(l),  ...,  i(m)  be  the  elements  of 
I  ,  in  increasing  order  (and  take  i(0)  »  0  ).  Thus  the  generating 
•  sequence '  „ 

«> >|.|« 

for  contains  (B*  )  terms.  On  the  other  hand,  this  sequence  is 

linearly  independent  since,  e.g., 

[0]D#()tt  -  el  . 

Note  that 

dim  (*k+1  e  «k)  - 

since  (a  €  z"  :  |a|  ■  k)  is  in  obvious  1-1  correspondence  with 

{e  e  s  |a|  <  k}  .  This  reaffirms  the  well  known  identity 

rm  +  k1  ,  £k  rm-l+r \ 
i  n  J  "r*0  '  m-1  '  * 

We  now  consider  the  LIPd^jT)  *■  LIP( *Jt,span( [t] )t6T)  with  T  a 
subset  of  .  We  call  T  correct  (for  *k  )  if  the  LIP(«k,T)  is 
correct.  Since  is  not  the  tensor  product  of  univariate  polynomial 

spaces,  it  seems  unlikely  that  we  could  employ  the  tensor  product  con¬ 
struct  to  obtain  correct  T's  .  Yet  it  is  possible,  as  the  following 
example,  two-dimensional  for  simplicity,  shows.  Recall  that  the  linear 
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projector  of  polynomial  interpolation  at  points  Uq, 
written  in  Newton  form  as 


'  uk 


with 

Therefore,  any  partial  sum 


£0  ♦uiIu0 . uiJ 


♦ui(x)  s*  (x-u0)...(x-ui_1>  ,  all  i  . 


can  be 


RI 


( i, j)€I 

of  its  tensor  product  with  Pv  for  some  point  sequence  vQ,  . ..,  v^ 

5 

is  also  a  linear  projector,  given  that  I  C.  {0,...,k}  .  The  range  of 

2 

Rj  is  somewhere  in  •  To  insure  that  it  is  actually  )  * 

choose  I  *  {(i,j)  :  i+j  <  k)  .  With  this  choice,  ran  Rj  C  *k  and 

equality  must  hold  since,  by  just  counting  terms,  we  see  that  ran  R. 

2  A 

has  dimension  (  k  )  which  is  dim  x^CR  )  .  It  is  now  a  nice  exercise 
to  verify  that,  for  this  choice  of  I 


Rj  solves  the  LIP(»k,T) 


with  T  s*  {(u^,v.)  €  R 

The  same  construction  works  in 


i+j<k}  . 

m  variables  and  so  provides  the 
only  standard  choice  of  correct  point  sets  T  for  x^r!0)  •  This  is 
the  siaplicial  choice  which,  up  to  a  linear  change  of  variables,  is 
T  -  t(u  ,v  ,...,wr)  €  H®  s  p+q+...+r<k] 


p'  qr 

(Wp) 

an  affine  change  of  variables 


with  (Up),  (Vp),  ...,  (Wp)  given  sequences  of  real  numbers.  Note  that 


x  >  Ax  +  b 

leaves  x^  invariant,  hence  leaves  invariant  the  collection  of  correct 
point  sets  T  for  x^  . 

More  general  correct  point  sets  can  be  generated  with  the  aid  of 
the  following  theorem  due  to  Chung  &  Yao  [1977],  -  and  here  I  must 
thank  Or.  A.  Genz  for  pointing  out  this  reference  to  me. 

Theorem  4.1.  If  the  point  set  T  c  rf**  has  cardinality 
dim  x^di?1)  ,  and,  for  every  t  €  T,  there  exist  k  distinct 
hyperplanes  on  which  all  points  in  T  lie  except  for  t  ,  then  T  is 
correct  for  x^d^1)  . 


It  is  clear  how  one  would  prove  this  theorem:  For  each  t  €  T  ,  we 
can  find,  by  assumption,  k  m-vectors  a j  ,  ...,  a^  and  scalars  bj  , 
. . . ,  bfc  so  that  the  k-th  degree  polynomial 

Ltix  | — >  (a1*x-b1)...(a)c*x-bJc) 

vanishes  on  T\t  but  not  at  t  .  This  implies  that,  for  any  given  g 
,  the  function 

£t€T  9(t>  VLt(t) 
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is  a  polynomial  of  degree  <  k  which  agrees  with  g  on  T  .  On  the 
other  hand,  the  'sequence*  ^t^teT  lineariy  independent  (since  it 
is  obviously  independent  over  T  ),  and,  by  assumption,  contains  exact¬ 
ly  | T |  »  dim  terms,  hence  must  be  a  basis  for  and  this 

establishes  the  uniqueness  of  the  interpolating  polynomial.  In  short, 
we  have  the  generalization  of  Lagrange's  way  of  treating  univariate 
polynomial  interpolation. 

A  particularly  striking  instance  of  such  a  correct  set  T  are  the 
"natural  lattices*  of  Chung  &  Yao  [1977)  rediscovered  recently  by 
Dahmen  &  Micchelli  [1980],  and  also  by  Hakopian  [1981] 2  •  Pick  n 
points  a^  ,  . . . ,  an  in  if  so  that  the  points  0  ,  a^  ,  . . . ,  an  are 
in  general  position,  i.e.,  any  m+1  of  them  are  affinely  independent. 
To  recall,  m+1  points  bQ  ,  ...»  bm  in  if  are  affinely  independent 
if  their  affine  hull  is  all  of  if  ,  i.e.,  if  volm  conv  (b^JJ  *  0  • 
Then,  for  any  subset  I  of  (l,...,n>  with  |I|  *  m  ,  there  exists 
exactly  one  Xj  for  which 

1  +  a^Xj  -  0  ,  all  i  6  I  . 

(Indeed,  since  0  ,  (ai)iei  are  affinely  independent,  the  sequence 
(a^)i£I  must  be  linearly  independent.)  Further,  for  this  Xj  ,  we  must 
have 

1  +  aj*Xj  /  0  ,  all 

(since  1  +  aj*Xj  »  0  implies  that  (ai>iej  »  ®j  *11  lie  in  the 
hyperplane  (x  6  if1:  1+x‘Xj  *  0}  ,  hence  are  not  affinely  independent 
which,  by  assumption,  is  possible  only  if  j  e  I  ).  We  conclude  that 
T  :»  [Xj  :  I  C  {l,...,n),  |l|  «  m} 
is  correct  for  w^df)  with  k  n-m  ,  since  |T|  »  (”)  *  dim  and 

L.:x  1 — >  1  f 

j*I 

with  Lj(Xj)  »  . 

It  is  a  nice  exercise  to  develop  a  Newton  form  for  the  resulting 
polynomial  interpolation  scheme.  This  leads  to  a  particular  generaliz¬ 
ation  of  divided  differences  quite  different  from  the  tensor  product 
construction  with  which  we  began  this  discussion. 

The  Newton  approach  has  recently  been  generalized  by  Gasca  & 

Maeztu  (1980).  Although  the  idea  is  proposed  in  if  ,  I  shall  follow 
its  authors  and  discuss  details  only  in  R2  .  Start  with  a  straight  line 
1  +  ai*x  *  0  .  (This  is  not  quite  the  most  general  line,  but  that 
doesn't  matter.)  Add  a  bunch  of  lines  1  +  a^*x  ■  0  ,  intersecting  the 
first  line  at  distinct  points  xlir  i>l .  Form  the  polynomials 
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plis*  *“ >  (1  +  all*x)  **•  {1  +  al,i-l*x)  '  • 

The  LlPtspantpjj^)  ,span(  [x^] ) )  is  correct  since  the  Gram  matrix 

(Pl j ( x i i > )  is  triangular  with  nonzero  diagonal,  hence  invertible. 

Now  add  a  second  line  1  ♦  *2**  “  ®  intersecting  the  first  line 
at  a  point  other  than  the  x^'s  (if  at  all),  and  add  a  second  bunch  of 
lines  1  +  a2^*x  *  0  ,  intersecting  the  second  line  at  distinct 
points  x2i,  i“l» • • • ,  m2  •  Form  the  corresponding  polynomials 

p2iSX  *“■ >  (1  +  ai*x)(1  +  a21*x)  (1  +  a2,i-l*x>  '  * 

Then  (and  this  is  the  salient  part  of  the  construction),  the  single 

linear  factor  (1  +  a2*x)  in  the  p2^  makes  them  vanish  at  all  the 
earlier  interpolation  points  Xjj  .  with  this,  the  matrix  (Prj(*si>) 
in  lexicographic  order  is  triangular  with  nonvanishing  diagonal,  and 
thus  the  LIP(span(prj ) ,span( [xsi] ) )  is  correct. 

The  general  pattern  is  now  clear.  What  is  less  clear  is  just 
what  span(pcj)  might  be  and,  in  more  than  two  variables,  things 
become  horrendous.  Still,  for  certain  regular  choices  (see  Maeztu 
[1982]),  span(prj)  can  be  shown  to  coincide  with  and  the  corre¬ 

sponding  correct  point  set  can  be  more  general  than  the  simplicial 
choice,  but  not  more  general  than  those  covered  by  Theorem  4.1. 


In  his  1978  thesis  (see  Kergin  [1978],  [1980]),  Paul  Kergin  pro¬ 
poses  a  totally  different  approach  to  multivariate  polynomial  interpol- 
ation  which,  in  a  way,  gave  impetus  to  all  the  material  yet  to  be  dis¬ 
cussed  in  these  lectures.  I  begin  with  Kergin’s  result  as  he  stated  it. 
Theorem  4.2.  For  any  point  sequence  t„  ,  . . . ,  t„  in  R®  there 


>int  sequence 


exists  exactly  one  map  PsC^^R®)  — >  x_  so  that 

— —  -■  **■  •  *  n  ■  ~  ■ 

( i)  P  is  linear; 

(ii)  V  g6C^n^  V  0<k<n  V  q^Sx^  homogeneous  of  degree  k 
V  JC  {0,...,n}  with  |J|  =  k+1  , 

qk(D)Pg  =  qJc(D)g 
at  some  point  in  conv(tj)j€J  . 

Here,  as  earlier,  conv  T  denotes  the  convex  hull  of  the  point 
set  T  .  Further,  qk  €  x^  is  homogeneous  of  degree  k  in  case 

qk  -  I  A  ()“  . 

|«|-k 

Consider  for  a  moment  the  special  case  m  ■  2  .  The  stated  requir¬ 
ement  for  k  »  0  forces  Pg  to  agree  with  g  at  each  of  the  t^'s  . 
For  k  *  1  ,  we  have 

qk(D)  *  A< 1 ,0 ) Dlg  +  A(0,l)D2g  ’ 
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Thus,  if  f  tj  ,  then  Pg  would  have  to  match  any  such  derivative 
q^(D)  somewhere  on  the  segment  between  t^  and  tj  .  This  condition 
is  already  satisfied  in  case  q^(D)  is  in  the  same  direction  as  the 
segment  since  Pg  matches  g  at  and  tj  .  Therefore,  this 

imposes  just  one  additional  condition,  viz.  that  the  derivative  normal 
to  the  segment  be  matched  somewhere  along  the  segment.  If  t^  ■  tj 
for  i  ^  j  ,  we  get  two  additional  conditions,  viz.  that  Pg  have  the 
same  tangent  plane  at  t^  as  does  g  .  In  other  words,  we  obtain  oscu- 
latory  interpolation. 

An  extreme  case  of  osculatory  interpolation  occurs  in  case  tg  - 
...  ■  tn  .  Now  Pg  is  necessarily  just  the  Taylor  expansion  of 
degree  <  n  for  g  at  tQ  . 

Kergin  begins  the  proof  of  his  theorem  with  the  observation  that 
P  is  necessarily  continuous  on  x  :■  C^n*(G)  for  any  bounded  G 
containing  tg  ,  ...,  tR  .  Since  P  is  linear,  this  requires  only  to 
show  that  p  is  bounded  (on  X  )  .  This  latter  fact  Kergin  shows  by 
observing  that,  by  assumption,  the  leading  coefficients  of  Pg  agree 
with  the  corresponding  normalized  derivatives  of  g  at  certain  points, 
hence  can  be  bounded  in  terms  of  lgix  .  He  then  considers  Pg  -  Lg  , 
with  Lg  the  leading  terms  of  Pg  just  estimated, and  observes  that 
the  leading  terms  of  the  resulting  polynomial  are  of  lower  order  and 
must  interpolate  to  the  corresponding  normalized  derivatives  of  g  - 
Lg  at  certain  points,  hence  can  be  bounded  in  terms  of  lg  -  Lgix  , 
therefore  in  terms  of  ■ g l x  .  The  inductive  argument  is  now  clear. 

In  consequence,  P  can  be  understood  entirely  from  its  action  on 
a  fundamental  subset  of  ,  i.e.,  a  subset  R  whose  finite  linear 

combinations  are  dense  in  .  Kergin  chooses  R  to  consist  of 

socalled  plane  waves  (F.  John)  or  ridge  functions  (C.  A.  Micchelli), 

R  »-  {gol  x  g  6  C(n)(R)  ,  X  €  (H?)'}  . 

Such  a  function  is  constant  in  all  planes  normal  to  a  certain 
direction.  Explicitly, 

(goA ) (x)  -  g(X*x)  -  g(l®X(i)x(i))  . 

Note  that  it  is  sufficient  to  take  just  one  suitable  g  ,  e.g., 
gst  I— >  eifc  . 

Next,  Kergin  shows  uniqueness.  To  be  sure,  the  claim  is  not  that, 
for  a  given  g  ,  there  is  a  unique  Pg  6  *n  satisfying  the  conditions 
described  in  (ii),  for  that  is  not  true.  For  example,  the  function 
g:x  | — >  x ( 1 ) x ( 2 )  has  all  functions  psx  | — >  ax(2)  with  0  <  a  <  1  as 


linear  "interpolants"  at  the  points  ( 0 « 0 )  and  (1>0)  in  that  sense. 
Rather,  Kergin  claims  the  uniqueness  of  such  a  linear  map  P  and 
proves  it  by  showing  that  plane  waves  (with  the  univariate  function 
g  a  polynomial)  have  unique  "interpolants*.  Given  the  many 
conditions  P  has  to  satisfy,  the  uniqueness  is  not  surprising.  The 
hard  part  is  to  show  existence. 

For  this,  Kergin  introduces  (in  rather  different  notation)  the 
linear  functionals 


l  sv 

/  g  :■  JQ  ...JQ  gtXQ+s^x^  ...  +  s^Vx^)  dsJc...ds1  (4.1) 

(*0 . *icl 

and  sets 

Q  :■  span  {  g| — >  ^  qk(D)g  :  |J|  »  k+1  ,  k*0,...,n)  . 

Here,  tj  :*  (tj)j€J  ,  with J J  C  (0,...,n)  as  before.  Then,  by  some 
hard  counting,  Kergin  shows  that  dim  Q  <  dim  »n  .  Add  to  this  the  fact 
that 

QiA*n  “  {0)  (4.2) 

and  you  can  conclude  that  the  LIP(*n,Q)  is  correct.  Now  take  for  P 
the  resulting  projector.  Then 

/[t  ]<J)c(D,(9  *  Pg)  ”  0  ' 

J 

hence  qk(D)(g  -  Pg )  -  0  at  some  point  in  conv  tj  .  The  claim  (4.2) 

is  established  by  an  inductive  arguments  If  p  €  »nA  Qi  ,  then,  for  all 

|a|  ■  n  ,  Jf  .  ,Dap  *  0  ,  therefore  Dap  *  0  ,  i.e.,  p  6w  ,  , 

n-j. 


Micchelli  &  Milman  (1980]  give  a  striking  formulation  of  Kergin’s 
interpolation  scheme  which  shows  it  to  be  a  "lifting"  of  the  Newton 
form  of  the  univariate  interpolating  polynomial.  Micchelli  came  to  this 
by  noticing  that  Kergin's  linear  functionals  (4.1)  are  closely  related 
to  the  divided  difference  (as  the  notation  used  in  (4.1)  already 
intimates)  via  the  Hermite-Genocchi  formula  (see  Nflrlund  (1923;  p.16]): 


Dkg 


tTo'”*'Tk^9  =  / 

lT0'*..' “k1 

for  any  sufficiently  smooth  univariate  g  ,  a  fact  easily  proved  by 
induction.  This  allows  us  to  write  the  univariate  polynomial  interpol- 
ant  PTf  in  Newton  form  as 


(Pt9)(x) 


Z0  (x"V 


(x“Tk-l) 


J  Dkg 

[XQ, . . . ,Tk) 


Also,  recall  that  [i-jpa  =  ( t  ] g  for  all  J  c.  (0,...,n)  .  Now 

J  T  J  — 

consider  the  Micchell i-Milman  definition 
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for  any  £ 
inuous  on 
set  K  of 

therefore 


Pf 


x  1 — > 


E  It  It]  DX-t  DX“t  ^ 

k*0  1  c0  ’  ‘ 4  * '  ckJ  x  c0  ^k-l 

€  X  s»  C(n)(RP)  •  The  resulting  map  P  is  linear  and  cont- 

X  ,  hence  can  be  understood  by  looking  at  its  action  on  the 

plane  waves.  For  f  =  gox  €  R  ,  one  computes 
,( 1) 


V 


I1?  y(i)  g'*'  (X**)X(i) 


(i*y) 


g(1>ox 


x-t. 


■'x-t 


and  so 


k-1 


X*(x- 


V 


X*(x-tk_1)  g(k)oX 


Pf(x)  -  l"  X*(x-trt)  ...  X*(x-tk  .)  /  g(k)  . 

00  K  1  U*t0,...,X*tk] 

The  last  integral  equals  [X*tQ , . . . ,X*tk)g  ,  by  the  Hermite-Genocchi 

formula.  Therefore ,  finally, 


P(goX)  -  (p(x*t.)g)oX  * 

This  is  the  crucial  observation.  It  shows  that  ran  PC.*  and  that, 

—  n 

for  any  f  »  gox  6  R  ,  any  polynomial  qk  homogeneous  of  degree  k 
and  any  JC  {0,...,n}  with  |J|  -  k+1  , 

/it.i  VD1M  *  <lt,i  V01'  • 

since,  for  such  an  f  , 


and 


qk(D)f 


I  AD1#,1)...  D“,n0  f 
|a|-k  °  1  m 


C  *.‘“9 


(k) 


oX 


^[tjl  ^(X*^)^ 


(k) 


oX 


|a|-k 

lX*VP(X*ti)g  “  * 


This  establishes  that  (4.3)  is  a  formula  for  Kergin's  map.  ||| 

Micchelli  [1980]  offers  additional  detail,  e.g.,  the  error  formula 
one  associates  with  the  Newton  form  which  leads  to  a  constructive  proof 
of  the  Bramble-Hilbert  lemma  (see  Bramble  &  Hubert  [VI  *01). 

Kergin's  scheme  raises  some  questions*  Ir.  contrast  to  its  univar¬ 
iate  antecedent,  it  requires  derivative  information  even  when  all  the 
interpolation  points  are  distinct.  This  has  been  remedied  recently  in 
Hakopian  [1981] ^  by  the  simple  device  of  lowering  the  order  of  the 
derivatives  appearing  in  the  definition  of  the  linear  functionals  which 
make  up  Kergin's  interpolation  conditions  Q  .  One  may  also  investigate 
other  variants  of  Kergin's  scheme.  In  particular,  Kergin's  scheme  makes 
it  very  attractive  to  study  the  linear  functionals  (4.1),  as  a  basis 
for  a  suitable  definition  of  the  divided  differences  of  a  function  of 
several  variables.  This  was,  in  fact,  the  consideration  which  led 
Kergin  to  his  scheme  in  the  first  place. 
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The  study  of  these  linear  functionals  led  C.  A.  Micchelli  to  the 
recurrence  relations  for  multivariate  B-splines  and  so  opened  up  that 
fruitful  area  of  research  which  is  the  topic  of  the  remaining  lectures. 


5.  Multivariate  B-splines 


Following  the  lead  of  Schoenberg  [1965) ,  the  multivariate  B-spline 
M( • | t0, . . . , tn)  was  defined  in  de  Boor  [1976]  by  the  rule 


M(  x  1 1. , . . . ,  t  ) 


voln-m^~lx^[t0>-;>tn1 
voin  ft0,...,tnJ 


€  R™ 


thus  generalizing  a  particular  characterization  of  the  univariate 
B-spline  due  to  Curry  &  Schoenberg  [1966).  Here,  t0,  . ..,  tn  are 
points  in  R°  ,  [K]  is  the  convex  hull  of  the  point  set  K  ,  and  P 

is  the  canonical  projector 

P  s  Rn — >  Rm  :  x  >  (x(i))®  . 

Further,  volk(K)  is  the  k-dimensional  volume  of  the  set  K  . 

Such  a  B-spline  is  a  nonnegative  piecewise  polynomial  function  of 
degree  «  k  :«  n  -  m  ,  its  support  is  [PtQ, . . . ,Ptn)  ,  and  it  is  in 
cn-m-l  ag  lQng  as  the  ptQ  ,  ptn  are  in  generai  position. 

All  this  will  be  proved  shortly. 

In  1978,  C.  A.  Micchelli  [1980]  proposed  an  equivalent  but  more 
suitable  and  flexible  definition  of  M( • | tQ , . . . , tn)  as  the 
distribution  c**  CgfR?*)  given  by  the  rule 


M(*|t-,...,t)  :  f  |— >  nl  / 


fop 


(5.1) 


[to'-  •  ‘'V 


This  definition  makes  sense  even  if  the  t^*s  are  not  in  general 
position.  M(* |t0,...,tn)  is  a  function  (in  L.(R®)  )  if  and  only  if 
volm(Pto, . . . ,Ptn)  *  0  and,  in  that  case, 

nl  /  fop  (5.2) 


"<-ito . V' 


lto"  *  •  'S*1 


More  than  that,  Micchelli  [1980)  proved  recurrence  relations  for 
these  multivariate  B-splines,  of  the  following  form. 


Theorem  5.1.  ( 1 )  If  x  =  I  °ipti  with 


*  0  ,  then 


DxM(,IV**-'tn)  =  n  1  *iM<’|t0 . ti_lfti+1 

( ii )  If  x  =  I  °ipti  with  la^  =  1  ,  then 

(n-m)  M(x|t0,...,tn)  >=  n  Z  a.M(  x  |  tQ, . . . ,  t  ,_1 ,  t.+1 


,...,tn) 
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These  recurrence  relations  were  proved  almost  simultaneously  by  a 
different  approach  by  Dahmen  (1979) ^  and  have  since  then  been  given 
different  proofs  by  Micchelli  (1979],  H&llig  (1980),  Hakopian  (1980), 
de  Boor  k  H&llig  (1981)  and  perhaps  others.  I'll  now  give  a  version  of 
this  last  proof,  from  de  Boor  t  H&llig  (1982),  not  only  because, 
naturally,  I  like  it  best,  but  because  it  covers  a  more  general 
situation  than  described  so  far. 

To  begin  with,  I  have  to  clear  up  an  inconsistency  in  the  notation 
I  have  employed.  This  inconsistency  shows  up  in  (S.2)  where  both 

!m  and  / 

®  (tjj,...,  tn) 

occur  and  (tQ,...,tn]  is,  off-hand,  not  meant  as  the  convex  hull  of 
the  points  tQ,  ...,  tR  ,  but  rather  as  an  indication  that  the  integral 
is  to  be  formed  as  described  in  (4.1).  These  two  meanings  only  differ 
by  a  scalar  factor,  though,  viz.  the  factor  nlvoln(tg, . . . ,tn)  ,  hence 
could  be  made  to  coincide  if,  in  (4.1),  we  multiplied  the  right  hand 
side  by  nlvoln  (tg,...,tn)  .  I  settle  this  inconsistency  instead  by 
abandoning  from  now  on  the  interpretation  (4.1)  and  entirely  rely  on 
the  more  naive  interpretation  of 

^(t0,...,tn) 

as  the  integral  over  the  convex  hull  of  the  points  tQ  ,  ...»  tn  . 

Consider  now,  more  generally,  a  polyhedral  convex  set  B  in  R11  , 
some  linear  map  p  into  r”1  and  having  B  in  its  domain,  and  the 
distribution  Mg  on  R®  defined  by  the  rule 

MB  :  f  Y—  >  JB  fop  ,  (5.3) 

i.e.,  as  the  P-shadow  of  B  .  The  simplex  spline  M( • | tfl, . . . ,tn) 
results  when  B  -  (tQ,.. ,tn)/voln(tg,.. . ,tn)  and  P  is  the  canonical 
map  Rn  — >  R™  ,  but  it  has  already  turned  out  to  be  fruitful  to  allow 
more  general  sets  B  and,  usually  as  a  simplification,  more  general 
linear  maps  P  . 

In  any  case,  the  recurrence  relations  can  be  proved  in  this 
generality.  The  relevant  observation  is  that  the  boundary  of  such  a 
convex  polyhedral  set  consists  again  of  convex  polyhedral  sets  B£  ,  of 
one  dimension  lower,  hence  Stokes'  Theorem  can  be  used  to  relate  M  : » 
Mg  to  Mi  !*  mb.  *  For  this,  we  also  need  the  corresponding  outward 
normal  to  B  at  B^  relative  to  the  affine  hull  of  B  ,  and  an 

arbitrary  point  bA  in  the  affine  hull  of  B^  .  For  simplicity,  we 
assume  that  B  is  a  body,  i.e.,  B  has  R*1  as  its  affine  hull.  With 
these  assumptions  and  notations,  the  following  theorem  holds. 
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Theorem  5.2.  (i)  DpzM  *  -  E^  z-v^  M^  ,  all  z  6  Rn  . 

(ii)  (n-m)M(Pz)  -  Ii  (bi-z)*vi  M^Pz)  ,  all  z  6  Rn  . 
(ill)  DM  -  (n-m)M  -  ti  M^  . 


The  proof  of  (i)  is  immediate: 


(DpzM)f  -  -  JB  (Dp2f)op  -  -  /B  Dz(fop)  -  -  fiB z*v  fop 

•  ~  zi  fB  **vi  f°P  “  “  ri  z*vi  Mif  ‘ 

This  uses  the  fact  that,  by  definition,  the  derivative  DyM  of  the 
distribution  M  is  the  distribution  obtained  by  the  rule 
f  I—  >  M(-Dyf)  ,  and  the  standard  interplay 

Dy(fop)  -  (Dpyf)op 

between  differentiation  and  linear  change  of  variables.  This  interplay 
also  proves  that 

( Df ) ( Px )  -  ( Dpx  f )  (  PX )  -  (Dx(fop))(x)  -  (d( fop ) ) ( x )  (5.4) 

with  D  the  first  order  differential  operator  given  by  the  rule 

r 

Df  S  F,  D,f 

r  3-1  1  1 

in  case  f  has  its  domain  in  R  ,  and  with 

( F j f ) ( x )  x( j )f (x)  ,  all  j  . 

Thus  (Df ) (x)  *  (Dxf)(x)  ,  and  the  adjoint  of  d  is  -E  d^F^  .  This  is 
of  use  in  proving  (iii):  We  have 


therefore 

and 


-  (dm) f 


DjFj  m  1+FjDj  , 

m 

-  /_(  E  D.F.fjop  -  mMf 

8  j-1  3  3 
n 

/_  E  D4F. (fop)  -  nMf 
8  i-1  1  1 


+  /B(Df)op 
+  fB  D(fop)  . 


Here,  the  last  integral  in  the  first  line  equals  the  last  integral  in 
the  second,  by  (5.4).  Therefore, 

n 

(DM) f  -  (n-m)Mf  -  E  /_  Dj F . ( fop)  . 

i-1  11 

This  settles  (iii)  since 

/b  DiFi(f0P)  -  if1  /8BV(i)  Fi(f0p)  “  ;3BC*v)  fop 

and,  on  the  facet  B^  ,  the  function  (»*v)  is  constant. 

Finally,  to  prove  (ii),  conclude  from  (i)  and  (iii)  that,  for 
any  z  with  Pz  *  x  , 
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0  -  (d  -  Dpz)M(X) 

-  (n-m)M(x)  -  I  ( x )  +  £  z^  Mi(x)  . 


As  an  exercise>  I  specialize  Theorem  5.2  to  the  situation  of 
Theorem  5.1.  This  means  that  B  »  [tg,...,tn]  and  that  we  may  set 
s *  [<tj)\ti]  ,  i-0,...,n  .  Then 

M  *  VOlnB  M(  »  |  tg , . .  .  ,  tR )  t 

therefore 


*  vo^n-l®i  *  I ^0 ' *  *  * ' ^i-1 '^i+1 ' *  *  * ,bn^  * 

Also, 

vol  B^  *  -  n  vol  B 

showing  that  the  coeffients  which  appear  in  (i)-(iii)  are,  for  this 
case,  closely  related  to  the  barycentric  or  areal  coordinates 
associated  with  the  simplex  B  .  In  any  case,  if  z  *  with 

I •  0  ,  then,  since  is  perpendicular  to  the  affine  hull  of 

(tj)  \t^  ,  we  have  ( t j— “  0  for  all  j  /  i  and  so 
z*v . 


Ij«j(tj-bi)*vi 


“i(ti‘bi>*vi  ' 


therefore,  from  (i), 

DPzM(,l(tj))  ”  DPzM/vo1  B  *  (-£i  Z*vi  Mj/vol  B 

■  n  ai  Mj/vol  Bj^  *  n  EiajM<  *  t  (tj)  \ tA ) 

This  proves  (i)  of  Micchelli's  Theorem  5.1,  with  x  *  Pz  .  For  Theorem 
5.1(ii),  we  have,  with  x  ■  Pz  ,  z  ■  and  Ej°j  “  1  »  that 


(bj-z)^ 


(bi”  Jajtj)*wi  “  Iaj<bi“tj)*vi 


«i(bi-t.)*vi  , 


therefore,  from  (ii), 

(n-m)M(x  |  ( tj ) )  -  £i(bi-z)*vi  M^xJ/vol  B 

*  (bi-ti)*vi  M^xJ/vol  B  -  £A  a^HUl  ( t j  )  Ntj)  . 


W.  Dahmen  has  pointed  out  to  me  that,  once  one  recognizes  that  the 
recurrence  relations  in  Theorem  5.1  can  be  written  as  in  Theorem  5.2  in 
terms  of  facet  normals,  then  Theorem  5.2  can  be  derived  from  Theorem 
5.1  using  the  fact  that  any  convex  polyhedral  body  can  be  triangulated, 
i.e.,  is  the  essentially  disjoint  union  of  simplices,  any  two  of  which 
are  either  disjoint  or  else  have  exactly  a  face  in  common. 

The  recurrence  relation  (iii)  was  first  stated  and  proved  (for 
simplices)  by  Hakopian  [1980],  in  his  simple  proof  of  Theorem  5.1. 


Particularly  useful  choices  for  B  other  than  the  simplex  include 
cones  and  boxes.  The  choice  of  the  cone 
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U0  +  Elai(ti~V  *  0  «  ^  .  all  i) 

with  vertex  tQ  and  generating  rays  t^-t0  ,  i-l,...,n  leads  to  the 
cone  spline  M  -  MB  .  This  is  the  truncated  power  introduced  and 
heavily  used  in  Dahmen  [1979],  in  direct  generalization  of  the 
functions 

R  — >  R  s  x  \—>  (x  -  t)+k_1 

familiar  from  univariate  spline  analysis.  Of  course,  every  such  cone 
spline  can  be  obtained  as  a  translate  of  the  p-shadow  of  the  standard 
cone 

R+n  <z  6  Rn  s  0  <  z(i)  ,  all  i  } 
for  an  appropriate  choice  of  the  linear  map  P  . 

Choice  of  the  box  or  parallelepiped 

{t0  +  E1  ai<ti_t0)  •*  0  *  «t  «  1  '  a11  *> 
gives  rise  to  the  box  spline  introduced  in  de  Boor  &  DeVore  [1981]  and 

further  studied  in  de  Boor  a  Hbllig  [1982] .  Any  such  box  spline  is  a 
translate  of  the  P-shadow  of  the  standard  box 

(z  €  Rn  :  0  <  z(i)  <  1  ,  all  i) 
for  an  appropriate  choice  of  the  linear  map  p  . 


Repeated  application  of  the  recurrence  relation  (i)  provides  the 
information  that,  for  arbitrary  vectors  y^  ,  yr  , 

D„  €  span{  M_  :  F  is  an  (n-r Jdim.face  of  B}. 

Yi  yr  B  F 

Thus  all  r-th  order  derivatives  of  Mg  are  in 


:  dim  PF  » 
MB  0  C(d) 


m 


V  (n-r) -dim. faces  F  of  B> 


L„  provided  PF  is 
m-dimensional  for  every  (n-r)-diraensional  face  F  of  B  .  This  allows 
the  following  conclusions,  asserted  earlier  for  the  simplex  spline: 
(iv)  MB  €  L^d) C  c(d-1)  ,  with 
d  :■  max  {  r 
and,  for  this  d  , 

-  M, 

<  k  :■  n-m  ,  with  supp  HbCFB  .  Indeed,  any  (n-k-l)-dimensional 
face  F  of  B  is  mapped  by  P  into  some  hyperplane  in  Rf11  ,  hence 
all  (k+l)st  order  derivatives  of 
some  hyperplanes. 
complement  of 


(v)  PB  _is^  m-dimensional,  then  Mg  is  a  pp  function  of  degree 

iny  (n-k-l)-di* 

R™ 

4g  have  their  support  entirely  on 
This  implies  that,  on  each  connected  component  of  the 


Mg  agrees  with  some  element  of 


{  PF  :  F  is  a  face  of  B 

. 


dim  PF  <  m  > 
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6.  Approximation  from  the  span  of  multivariate  B-splines 


A  solitary  B-spline  is  of  little  use  in  approximation.  Thus  we 
consider  now  a  whole  collection  B  :■  (B)  of  polyhedral  convex  bodies 

in  R0  ,  each  giving  rise  to  its  P-shadow  mb  ,  and  ask  just  how  B 

should  be  chosen  so  that  we  get  a  "useful"  family  (Mg)  .  We  use  the 
well  known  properties  of  Schoenberg's  univariate  B-spline  (see,  e.g., 
de  Boor  [1976])  as  a  guide. 

A  first  useful  property  of  univariate  B-splines  is  that,  properly 
normalized,  they  form  a  partition  of  unity.  This  is  not  so  hard  to 
achieve  with  multivariate  B-splines.  We  take  for  p  the  canonical 
map  R°  —  >  R®  and  choose  B  as  a  partition  of  some  slab  Rf®*A  in 
B11,  for  some  A  C  R*  (with  k  t>  n-m  ,  as  before).  Then 

t  Mg(x)  •  t  volj(.p“^x  AB  ■  volJtP“1xnljB  ■  vol^A  , 

B  B  B 

i.e.,  we  have  a  partition  of  the  constant  vol^A  •  Choosing  A  to  have 

volj^A  »  1  or  else  dividing  each  Mg  by  vol^A  ,  gives  the  desired 

partition  of  unity.  We  conclude  at  once  that  a  continuous  function  f 
can  be  approximated  from 


to  within 


SB  i"  span  (Mg)ggg 


with 


e(f l | B | )  , 


l»l  *•  •UPBBB  diaB  PB  • 

The  simple  approximation 

l  f(Tg)  Mg 

with  tq  €  B  ,  all  B  ,  is  that  accurate  (exercise). 

Can  we  do  better  in  case  f  is  smoother?  From  the  univariate 
theory,  we  would  expect  to  get 

dist( f ,SB)  -  0(|B|*+1)  (6.1) 

in  case  f  €  .  This  we  could  conclude  at  once  if  we  had  available 

a  quasi-interpolant  q  for  SB  ,  i.e.,  a  bounded  linear  map  Q  into 
SB  which  is  local  and  reproduces  .  A  typical  specification  of 
■local"  would  be  to  require  that  Qf|C  depend  only  on  f|N(C)  with 

N(C)  «•  U(  PB  s  PB  fiC  *  0  )  . 

Then,  as  in  Lebesgue's  inequality,  we  could  conclude  that 
f  -  Qf  -  f-p  -  Q(f-p)  ,  for  all  p  e  ix  , 

therefore 


Mf-Qf) |C«  <  l(f-p)|C»  *  lOll(f-p))N(c)i  , 
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and  so 

l(f  -  Qf ) |CI  <  (1  +  IQl)  distN(C)(f,*k)  . 

This  leads  to  (6.1). 

A  first  requirement  for  such  an  argument  is  that  c.  SB  .  This 
was  established  by  Dahmen  in  (197912  by  a  clever  argument  for  the  case 
that  B  is  a  triangulation  of  such  a  slab 

s  »-  bP*  a  , 

i.e.,  B  consists  of  simplices.  In  effect,  he  deforms  the  slab 
appropriately  to  S'  so  that  the  function 

p  j  x  h->  voljj  P-1x  hS' 

is  (locally)  a  polynomial  of  degree  k  ,  yet  S'  is  still  triangulated 
by  the  simplices  B’  ,  with  B*  the  simplex  into  whose  vertices  the 
deformation  sent  the  vertices  of  B  .  The  deformation  only  takes  place 
in  directions  perpendicular  to  sP  .  Therefore  p  carries  the  vertices 
of  B  and  the  corresponding  vertices  of  B'  to  the  same  points. 
Consequently, 

MB/volnB  ■  MB,/volnB'  , 

and  so 

p  -  EB6B(vol  B'/vol  B>  Mb  . 

He  is  able  to  modulate  the  deformation  sufficiently  to  obtain  a  basis 
for  in  this  way. 

This  description  neatly  avoids  discussion  of  some  very  nontrivial 
details. 

In  (1980]2  •  Dahmen  uses  A  *  [0,11*  .  There,  he  uses  his  result 
to  support  the  claim  that  (6.1)  holds  (called  Theorem  3.1  there),  even 
without  benefit  of  a  quasi-interpolant,  but  I  cannot  follow  the 
argument  (even  though  it  is  a  generalization  of  a  univariate  argument  I 
once  made  up).  I  must  therefore  go  the  way  of  quasi-interpolants.  Here 
I  come  across  a  difficulty  first  observed  by  Dahmen  (1980) ,  :  With  the 
choice  A  x  [0,1]  ,  the  B-splines  (Mg)  are  linearly  dependent. 

Dahmen  overcomes  this  difficulty  through  an  averaging  process  which 
provides  a  quasi-interpolant  bounded  in  terms  of  a  local  mesh  ratio. 

There  is  another,  practical  difficulty,  though,  with  this  setup. 
The  B-splines  MB  ,  their  support,  their  sets  of  discontinuities,  all 
depend  on  the  triangulation  B  of  the  slab  S  .  In  practice,  one  would 
like  to  start  with  some  triangulation  T  of  uP  and  then  construct 
B-splines  in 

*c  i-  »  AC(r) 

*k,T  *  k,T  C  * 

the  space  of  all  pp  functions  of  degree  <  k  and  smoothness  r 
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associated  with  that  triangulation.  Further,  one  would  like  to  find 
enough  B-splines  to  staff  a  basis  for  iF  _  ,  just  as  one  is  able  to  do 
for  m  ■  1  (which  acounts  for  the  "B"  in  their  name). 

This,  as  it  turns  out,  is  too  much  to  ask  for,  for  various 
reasons.  At  present,  the  nature  of  the  spaces  «F  _  is  not  at  all 
understood,  except  when  r  <  1  .  We  don't  even  know  exactly  when  we  are 
entitled  to  a  local  basis  for  that  space.  We  must  therefore  be 
satisfied  merely  to  construct  M.  6  iF  .  .  But  even  that  is  not 

°  JC  t  T 

possible  in  general,  since  the  sets  of  discontinuities  of  Mg  are 
associated  with  the  projections  of  faces  of  B  ,  hence  not  arbitrarily 
choosable.  In  particular,  even  if  B  is  a  simplex,  the  connected 
components  of  the  complement  of  the  discontinuity  set  need  not  be 
simplices. 

One  can  hope,  though,  to  obtain,  for  a  given  triangulation  T  ,  a 
linearly  independent  collection  (Mg)  of  B-splines  of  prescribed 
smoothness,  whose  subdivision  is  a  refinement  of  T  and  whose  span 
Sg  contains  .  It  is  not  difficult  to  do  this  for  the  univariate 
B-splines.  Take  for  A  the  standard  simplex  [e0,...,ek]  in  R* 

(with  e0  s*  0  ,  ej  :=  (4^)  for  j  >  0  )  and  triangulate  R*A  by  the 
simplices 

xi  ti+k*eli+k] ^  »  all  i  , 

with  (ti)  the  given  knot  sequence  in  R  and  [j)  the  remainder  on 
division  of  j  by  k+1  .  it  is  worthwhile  to  visualize  this 
construction  for  k  ■  1  and  2  .  The  construction  uses  the  total 
ordering  of  R  in  an  essential  way,  hence  it  is  not  easily  generalized 
to  the  general  case.  Nevertheless,  a  construction  of  the  desired  type 
was  found  by  Dahmen  &  Micchelli  11982)  and  independently  by  HOllig 
(1982).  The  underlying  idea  is  to  get  a  B-spline  basis  (Mg)  for 
Wj^T  and  then  to  "pull  apart  their  knots*.  Explicitly,  with 

A  -  (e0,...,ek) 

the  standard  simplex  in  R*1  and 

T  8“  lv0 . vmJ 

any  particular  simplex  in  the  triangulation  T  ,  the  simploid 

TxA 

is  triangulated  using  the  combinatorial  product  toa  .  I  find  it  helpful 
to  visualize  this  construction  in  the  following  way: 

In  the  cartesian  product  (vi)*(e^)  ,  there  are  (mjjj*t) 
nondecreasing  "paths"  with  endpoints  <*0,eQ)  and  (Tm,ek)  .  a  typical 
"path"  is  shown  in  the  following  figure. 
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Associate  with  each  such  "path*  the  simplex  c  in  t*A  spanned  by  the 
points  on  that  "path".  It  is  a  nice  and  worthwhile  exercise  to  show 
that  the  resulting  collection 

IT  -  <«) 

of  simplices  forms  a  triangulation  for  t*A  .  More  than  that,  if  we 
carry  out  this  construction  for  each  t  €  T  ,  making  certain  only  that 
the  ordering  of  the  vertices  of  each  t  used  is  consistent,  i.e., 
comes  from  a  total  ordering  of  all  the  vertices  appearing  in  T  ,  then 
we  obtain  a  triangulation 

r  t-  TOA 

for  RmxA  .  This  is  straightforward  except,  perhaps,  for  the  assertion 
that 

oAo'  -  [V9nv#l]  ,  all  o,  o'  6  E  , 

with  V  the  set  of  vertices  of  the  simplex  o  . 
o 

Now  consider  (M0)oeJ.  .  For  each  a  €  Ex  ,  Mo  has  support  in  t  , 
and  PV0  -  VT  ,  hence  M0  has  t  as  its  support  and  agrees  with  some 
polynomial  there.  For  a  given  t  6  T  ,  there  are  exactly  t1"**)  such 
Mff  ,  i.e.,  exactly  enough  to  staff  a  basis  for  .  since,  by 

Dahmen's  result,  c  Sj.  ,  it  follows  that  these  Hg  form  a  basis  for 
wk|t  *  We  conc^ude  that  (M0)0g£  is  a  basis  for  T  . 

It's  time  to  pull  apart  the  knots.  For  this,  denote  by 
V(  I  )  ■  V(T)  X  {fig  f  •  •  • 

the  vertices  of  the  triangulation  E  »  toa  .  Let 

F  i  V{ E )  — >  Rns  (v.ej^)  >  (Fv,ei)  . 

The  only  requirement  HSllig  [1982]  makes  on  F  (or  F  )  is  that  it  be 
locally  finite,  i.e.,  |ran  F  A  C|  <  •  for  all  bounded  sets  C  . 

Dahmen  &  Micchelli  [1982]  require  that  sign  det  o  ■  sign  det  Fa  , 
all  o  e  E  .  Here, 

det  e 

is  the  signed  volume  of  the  simplex  o  .  its  signature  depends,  of 
course,  on  the  ordering  of  its  vertices  and  the  ordering  is  meant  to  be 
the  ordering  in  V(E)  . 
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With  this,  I  am  ready  to  state  the  basic  result  of  this 
construction  due  to  Dahmen  a  Micchelli  and  to  HOllig. 


Theorem  6.1«  For  all  y  €  Rf"  , 

(1  +  x-y)*  -  I  c  (y)  M_  (x) 

o€I 

for  all  x  at  which  all  the  simplex  splines  occuring  on  the  right  are 
continuous,  with 

Cg(y)  s*  fci  sign(det  a)  det  GyFo 

and 

GyjR1”*^  — >  Rn:  <x,u)  |— >  (x,  (l+x*y)u)  . 


In  effect,  Gy  carries  out  Dahmen' s  appropriate  deformation  of 
the  slab  sP*xa  mentioned  earlier. 

Since  y  is  arbitrary,  we  conclude  that 

*k  £  SPI  * 

The  quasi-interpolant  is  now  immediate,  in  case  F  is  not  too 
violent.  Specifically,  assume  that,  for  each  t  €  T  ,  there  is  some 
ball  bT  which  is  contained  in  the  support  of  every  Mpa  with 
c  €  It  and  is  outside  the  support  of  every  other  .  This  implies 

that 


SpE  -  span  (Mpo)s6£^  on  bT 

and,  as  there  are  just  enough  simplex  splines  Mpo  to  staff  a  basis 
for  ,  this  implies  that 


mf«  |bx  8  %  ibT  for  every  •  «  *t  • 

This  allows  the  construction  of  linear  functionals  as 

normpreserving  extensions,  to  L^(bt)  say,  of  the  coordinate 
functionals  u0  on  which  carry  the  polynomial  MFojb  to  ®oai 

o'  e  I  .  Since  T 


all 


this  shows  that 


■upp  xa  -  bT  '  a11  0  8  JT  • 

(*0)  is  dual  to  (MPJ  ,  i.e.. 


Vfo- 


Fe ; 

aoo*  •  a11  081 


The  resulting  quasi-interpolant 

0  ,  f  V->  I  Xat  Hrg 

is  therefore  even  a  linear  projector  onto  Sp£  ,  i.e.,  it  reproduces 
all  of  Sp£  and  not  just  .  The  only  concern  is  its  size.  IQI  can 
be  bounded  in  terms  of  the  relative  size  of  b  in  t  .  in  particular. 


-34 


all  is  well  in  case  the  vertex  perturbation  map  F  is  not  too  violent. 

Theorem  6.1  provides  a  generalization  of  Marsden's  identity.  It  is 
fair  to  say  that  it  is  based  on  a  two-dimensional  version  of  Marsden's 
identity  first  proved  by  Goodman  &  Lee  [1981].  These  authors  provide, 

in  particular,  the  more  explicit  formula 

n 

C  (y)  «  kl  TT  (1  +  zs*y>  * 

j=3  3 

with  the  Zj  determined  as  follows:  The  vertices  Vj  of  Fo  are  of 
the  form 


PV.  x  ei 

i  *  0,...,k  »  n-2  ,  there  is  at  least  one 
.  This  leaves  exactly  two  possibilities: 


j  so  that 


and,  for  each 
v j  *  (PVj,ei) 

(i)  for  some  i  ,  there  are  three  vertices  of  the  form  (PVj,e^) 
Then  we  will  call  them  vQ,  v^,  v2  . 

(ii)  for  two  values  of  i  ,  there  are  two  vertices  of  the  form 
(PVj,e^).  Then  we  will  call  one  pair 
With  this,  we  take 

Pv  j 


'j 


‘j 


Pv, 


aff  [PvQ,Pv^]  Aaff  [Pv2,Pv3J 


,  if  (i)  | 

,  if  (ii)/' 


j  >  3 


j 


HOllig  [1982]  gives  a  simple  example  to  show  that  such  a  nice 
formula  with  linearly  factored  coefficients  is,  in  general,  not  to  be 
expected  for  m  >  2  .  Still,  for  the  practically  important  case 
m  ■  2  ,  these  simple  formulae  lead  Goodman  &  Lee  to  the  intriguing 
generalization 

V  :f  I — >  Z  f(t0)  MFo 

of  Schoenberg's  variation  diminishing  spline  operator.  This  operator 
V  is  obviously  positive  regardless  of  the  choice  of  the  t0  .  Goodman 
4  Lee  choose 

t0  (z3  +  ...  +  zn)/(n-2) 

since,  in  light  of  Theorem  6.1,  this  implies  that  V  reproduces  . 
They  are  able  to  prove  that,  for  any  continuous  f  ,  vf  converges  to 
f  in  the  uniform  norm  as  |FZ|  — >  0  ,  provided  only  that  all  the 
Zj's  for  o  lie  in  PFc  ,  an  issue  only  for  z3  and  only  in  case 
(ii). 


There 
of  certain 


is  an  analogous  quasi-interpolant  construction  for  the  span 
translates  of  a  box  spline  in  de  Boor  &  HOllig  [1982].  The 
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arguments  have  a  different  flavor,  though,  since  the  resulting  pp 
functions  have  regular  meshes,  hence  are  amenable  to  "cardinal  spline 
techniques  familiar  from  Schoenberg  [1973]. 


7.  Epilogue 

In  these  lectures,  I  have  touched  on  only  very  few  questions  of 
current  interest  in  multivariate  approximation  theory.  Even  if  I 
restrict  attention  to  splines  and  pp  functions,  there  are  several  areas 
of  current  research  which  I  had  intended  to  discuss  when  I  first 
prepared  for  these  lectures  but  which,  in  the  end,  I  did  not  manage  to 
fit  into  the  allotted  time. 

The  nature  of  the  space 

_  r 

S  s *  *.  „ 

k,T 

of  smooth  pp  functions  on  a  given  partition  T  of  some  Q  O  R  is  not 
at  all  understood.  Questions  of  interest  concern  the  existence  of  a 
locally  supported  basis  for  S  ,  the  dimension  of  S  ,  the  dimension  of 
the  subspace  consisting  of  those  f  6  S  which  vanish  (r+l)-fold  at  the 
boundary  of  Q  ,  the  degree  of  approximation  achievable  from  S  .  There 
is  the  conjecture  that  sufficiently  smooth  functions  can  be 
approximated  to  within  0(|T|S)  if  and  only  if  S  contains  a  local 
partition  of  every  p  €  «s_^  ,  but  attempts  to  prove  this  by 
construction  of  a  quasi-interpolant  have  required,  in  addition,  some 
kind  of  stability  of  the  partitions.  Work  concerning  dim  s  has  been 
done  only  for  m  -  2  and  initially  only  for  r  -  1  ,  the  first 
nontrivial  case.  See  Strang  [1974],  Morgan  &  Scott  [1975],  and  the 
survey  of  Schumaker  [1979].  Most  recently,  Chui  &  Wang  [1981] 3  have 
given  precise  results  for  certain  T  and  arbitrary  r  .  The  existence 
of  a  local  partition  of  unity  in  S  is  taken  up  in  de  Boor  6  DeVore 
[1981]  for  certain  regular  T  in  order  to  understand  better  the  degree 
of  approximation  from  S  .  These  questions  are  further  pursued  in  de 
Boor  a  Hflllig  [1982].  In  both  papers,  the  relationship  between  S  and 
B-splines  in  S  is  explored,  but  even  in  the  context  of  a  simple  and 
regular  T  (e.g.,  a  rectangular  grid  with  all  north-east  diagonals 
drawn  in),  this  relationship  is  not  yet  fully  understood. 

The  adaptive  choice  of  the  partition  T  is  the  topic  of  de  Boor  & 
Rice  [1979].  Dahmen  [1982]  describes  one  way  to  use  simplex  splines 
adaptively.  The  degree  of  approximation  achievable  from  S  by  proper 
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choice  of  T  is  the  topic  of  Dahmen,  DeVore  &  Scherer  [1980). 

r  2 

Practial  aspects  of  approximation  from  S  =  T  on  R  are  the 
topic  of  the  two  survey  papers  Barnhill  [1977)  and  Schumaker  [1976).  An 
interesting  comparison  of  methods  is  given  in  Franke  [1982).  And  then 
there  is  the  vast  literature  on  the  constructive  aspects  of  the  Pinite 
Element  Methodl  Some  references  of  particular  interest  to  Approximation 
Theory  are  :  Ciarlet  &  Raviart  [1972]^  2»  Courant  [1943),  Fix  &  Strang 
[1969),  Guglielmo  (1969)  ,  Strang  &  Fix  [1973],  but  this  is  clearly 
just  a  taste. 

The  variational  approach  to  splines  is,  of  course,  not  restricted 
to  the  univariate  situation.  Already  Golomb  &  Weinberger  [1959) 
consider  particular  bivariate  examples  as  illustrations  of  the  general 
theory.  This  theory  has  the  following  setting.  A  collection  A  of 
continuous  linear  functionals  on  some  linear  space  X  is  given.  Since 
the  problem  will  only  involve 

A  :■  {x  e  X  :  Xx  *  0  for  all  IGA)  , 
we  might  as  well  assume  that  A  is  a  closed  subspace  of  X  .  Further, 
a  bounded  linear  map  T  from  X  to  some  normed  linear  space  Y  is 
given.  The  problem  is  to  determine,  for  given  x  6  X  ,  if  possible,  an 
element  x*  at  which  the  map 

x  +  Ax  — >  R+s  y  | — >  iTyi 

takes  on  its  minimum.  Such  a  minimizer  x*  is  called  a  (T,A)-spline 
interpolant  to  x  .  The  word  " interpolant"  is  appropriate  since  x* 
agrees  with  x  on  A  .  Golomb  6  Weinberger  [19591  deal  with  the 
special  case:  X  ■  Y  is  a  Hilbert  space  and  T  =  1  .  In  this  setting, 
the  map  x  | — >  x*  is  just  the  orthogonal  projector  onto  A  considered 
as  a  subspace  of  X  .  This  schizophrenic  nature  of  the  interpolation 
conditions  A  ,  linear  functionals  on  the  one  hand  and  elements  of  X 
on  the  other,  is  at  the  heart  of  the  practical  application  of  this 
theory.  In  standard  Hilbert  spaces  of  smooth  functions  on  some  domain 
0  ,  the  linear  functional  [t)  of  evaluation  at  some  point  t  turns 
into  the  section  G(t,*)  of  the  appropriate  Green’s  function.  In 
particular,  when  X  =  L2^^[a,b]  ,  then  [t]  is  represented  by  a  pp 
function  of  degree  2k-l  and  in  c<2k-2)  with  just  one  breakpoint, 
at  t  . 

I  was  held  back  from  exploring  multivariate  splines  obtained  in 
this  way  by  the  realization  that  this  would  require  me  to  obtain  and 
work  with  the  Green's  function  relevant  to  X  .  This  would  usually  not 
be  polynomial  nor  locally  simple,  and  would  depend  essentially  on  the 


domain  Q  .  Duchon  [1976],  [1977]  dealt  with  such  objections  by  the 
very  effective  device  of  choosing  all  of  R™  for  Q  .  The  resulting 
thin  plate  splines  have  already  found  practical  use.  Their  theory  is 
described  invitingly  in  Meinguet  [19791. 
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